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Abstract 

The small wavenumber k behavior of the structure factor S(k) of overcompressed amorphous 
hard-sphere configurations was previously studied for a wide range of densities up to the maximally 
random jammed state, which can be viewed as a prototypical glassy state [A. Hopkins, F. H. 
Stillinger and S. Torquato, Phys. Rev. E, 86, 021505 (2012)]. It was found that a precursor to 
the glassy jammed state was evident long before the jamming density was reached as measured by 
a growing nonequilibrium length scale extracted from the volume integral of the direct correlation 
function c(r), which becomes long-ranged as the critical jammed state is reached. The present 
study extends that work by investigating via computer simulations two different atomic models: 
the single-component Z2 Dzugutov potential in three dimensions and the binary-mixture Kob- 
Andersen potential in two dimensions. Consistent with the aforementioned hard-sphere study, we 
demonstrate that for both models a signature of the glass transition is apparent well before the 
transition temperature is reached as measured by the length scale determined from from the volume 
integral of the direct correlation function in the single-component case and a generalized direct 
correlation function in the binary-mixture case. The latter quantity is obtained from a generalized 
Orstein-Zernike integral equation for a certain decoration of the atomic point configuration. We 
also show that these growing length scales, which are a consequence of the long-range nature of 
the direct correlation functions, are intrinsically nonequilibrium in nature as determined by an 
index X that is a measure of deviation from thermal equilibrium. It is also demonstrated that 
this nonequilibrium index, which increases upon supercooling, is correlated with a characteristic 
relaxation time scale. 
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I. INTRODUCTION 



A quantitative understanding of nature of the physics of the glass transition is one of 
the most fascinating and challenging problems in materials science and condensed-matter 
physics. A sufficiently rapid quench of a liquid from above its freezing temperature into a 
supercooled regime can avoid crystal nucleation to produce a glass with a relaxation time 
that is much larger than experimental time scales, resulting in an amorphous characteristic 
state (without long-range order) that is simultaneously rigid 1 . A question that has received 
considerable attention in recent years is whether the growing relaxation times under super- 
cooling have accompanying growing structural length scales. Two distinct schools of thought 
have emerged to address this question. One asserts that static structure of a glass, as mea- 
sured by pair correlations, is indistinguishable from that of the corresponding liquid. Thus, 
since there is no signature of increasing static correlation length scales accompanying the 
glass transition, it identifies growing dynamical length scales.- - - The other camp contends 
that there is a static growing length scale of thermodynamic origin^ and therefore one need 
not look for growing length scales associated with the dynamics. 

In the present paper, we employ both theoretical and computational methods to study 
two different atomic glass-forming liquid models that support an alternative view, namely, 
the existence of a growing static length scale as the temperature of the supercooled liquid 
is decreased that is intrinsically nonequilibrium in nature. This investigation extends re- 
cent previous work 7 in which this conclusion was first reached by examining overcompressed 
hard-sphere liquids up to the maximally random jammed (MRJ) state.- (For a hard-sphere 
system, compression qualitatively plays the same role as decreasing the temperatnre in an 
atomic or molecular system; see Ref. |9|.) The MRJ state under the strict-jamming con- 
straint is a prototypical glass in that it lacks any long-range order but is perfectly rigid such 
that the elastic moduli are unbounded.- 1 ^ This endows such packings with the special hy- 
peruniformity attribute. A statistically homogeneous and isotropic single-component point 
configuration at number density p is hyperuniform if its structure factor 

S(k) = l + ph(k) (1) 

tends to zero as the wavenumber k — > 0,— where h(r) = Q2{r) — 1 is the total correlation 
function, g2{r) is the pair correlation function, and h(k) is the Fourier transform of h(r). 
This condition implies that infinite-wavelength density fluctuations vanish. 



It was theoretically established that hyperuniform point distributions are at an "inverted" 
critical point in that the direct correlation function c(r), rather than the total correlation 
function h(r), becomes long-ranged, i.e., it decays more slowly than —1/r in <i-dimensional 
Euclidean space M. d , where r is the radial distance.— The Fourier transform of direct corre- 
lation function c(k) is defined via the Ornstein-Zernike equation: 12 

c(k) ~ ~ Kk) ~ m - 1 (2) 
{ ) ~ 1 + ph(k) " pS(k) ■ (2) 

It is immediately clear from this definition that the real-space volume integral of the direct 
correlation function c(r) diverges to minus infinity for any hyperuniform system, since the 
denominator of © tends to zero, i.e., 

limc(fc) = / c(r)dr — > — oo (3) 
fc^o J Rd 

MRJ packings of identical spheres possess a special type of hyperuniformity such that S(k) 
tends to zero linearly in k as k — > 0, implying quasi-long-ranged negative pair correlations 
(anticorrelations) in which h(r) decays as a power law —1/r 4 or, equivalently, a direct 
correlation function that decays as c(r) ~ —1/r 2 for large r, as dictated by Eq. (2).— These 
anticorrelations reflect an unusual spatial patterning of regions of lower and higher local 
particle densities relative to the system density. This quasi- long-range behavior of h(r) is 
distinctly different from typical liquids in equilibrium, which tend to exhibit more rapidly 
decaying pair correlations, including exponential decays. 

Reference |7| examined overcompressed hard-sphere configurations that follow Newtonian 
dynamics for a wide range of densities up to the MRJ state. A central result of this study 
was to establish that a precursor to the glassy jammed state was evident long before the 
MRJ density was reached as measured by an associated growing length scale, extracted 
from the volume integral of the direct correlation function c(r), which of course diverges 
at the "critical" hyperuniform MRJ state. It was also shown that the nonequilibrium sig- 
nature of the aforementioned quasi-long-range anticorrelations, which was quantified via a 
nonequilibrium index X, emerges well before the jammed state was reached. 

These results for nonequilibrium amorphous hard-sphere packings suggest that the direct 
correlation function of supercooled atomic models in which the atoms possess both repulsive 
and attractive interactions should provide a robust nonequilibrium static growing length scale 
as the temperature is decreased to the glass transition and below. Here we show that this is 



indeed the case by extracting length scales associated with standard and generalized direct 
correlation functions. In particular, we study the single- component Z2 Dzugutov potential 
in three dimensions and the binary-mixture Kob- Andersen potential in two dimensions. The 
Z2 Dzugutov potential for a single-component many-particle system in three dimensions has 
the following form:— 



The first term in (J3J) models Friedel oscillations for a metal with Fermi wave vectors of mag- 
nitude hp, while the second term adds a strong repulsion for sufficiently small interparticle 
separations. The parameters a and b control the relative strengths of both contributions and 
define the energy scale. The cutoff r c is selected to be at the third minimum of the potential, 
while the constant Vq is present to make the potential continuous at the cutoff. The param- 
eters 77, a, and n control the shapes of both functions in (j4]). The Kob- Andersen model for 
a two-dimensional binary mixture is given by a truncated Lennard- Jones potential: 15 



The parameter e aj 3 controls the strength of the attraction between two particles of species 
a and j3, while a a p is equal to 2" 1 / 6 times the distance between both particles at which the 
attraction is maximal. 

It is known that overcompressing a hard-sphere system is analogous to supercooling a 
thermal liquid, but to what extent does this analogy hold? Roughly speaking, a rapid 
densification of a monodisperse hard-sphere system leads to the terminal MRJ state (with 
packing fraction of about 0.64), which we have noted is a prototypical glass.- At this singular 
state, the system is never able to relax and hence the associated relaxation time is infinite.— 
Slower densification rates lead to other jammed states with packing fractions higher than 
0.64. 9 ' Moreover, it has been shown that below 0.64, metastable hard-sphere systems have 
bounded characteristic relaxation times,— ^ including the range of packing fractions of about 
0.58 ~ 0.60 (depending on the densification rate) that has been interpreted to be the onset 
of a kinetic glass transition.-^ Above a particular hard-sphere glass-transition density, the 
system is able to support a shear stress on time scales small compared to a characteristic 




(4) 
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relaxation time. Clearly, increasing the density of a hard-sphere system plays the same role 
as decreasing temperature of a thermal liquid. In a thermal system, a glass at absolute zero 
temperature has an infinite relaxation time classically, and hence this state is the analog 
of the hard-sphere MRJ state. The glass transition temperature T g , which depends on the 
quenching rate and possesses a bounded characteristic relaxation time, is analogous to the 
aforementioned kinetic transition in hard spheres. These strong analogies between glassy 
hard-sphere states and glassy atomic systems lead one to believe that the results of Ref. ?] for 
the former extend to the latter. Indeed, here we demonstrate that the aforementioned length 
scales grow as as the temperature is decreased to the glass transition T g and below. Moreover, 
we show that the nonequilibrium index X, previously shown^ to increase as a hard-sphere 
system is densified to the MRJ state, also grows for T < 2.2T g . This nonequilibrium index 
is also shown to be correlated with an early relaxation time r. 

In Sec. [Til we introduce two generalizations of the direct correlation function c(r) which 
apply for two-component systems. In Sec. II I II we describe the numerical techniques and 
parameters used in our simulations, while in Sec. IIVI we present the results we extract 
from these simulations. The latter includes the demonstration of the existence of growing 
nonequilibrium static length scales upon supercooling the two atomic-liquid models that we 
consider. Moreover, we show that the nonequilibrium index X is positively correlated with 
an early relaxation time, both of which increase as the temperature is decreased to the glass 
transition temperature and below. We conclude in Sec. |V] with a summary of our results 
and of their impact. 



II. STRUCTURAL SIGNATURES OF LARGE- WAVELENGTH DENSITY FLUC- 
TUATIONS IN BINARY MIXTURES 

It has been shown that for maximally random jammed binary sphere packings, the 
standard structure factor S(k), determined from the particle centroids, cannot be used 
to ascertain whether the system is hyperuniform, unlike the single-component MRJ sphere 
packing.— 1 ^ Instead it was shown that the spectral density defined below, can be 

employed to determine whether a binary MRJ packing is hyperuniform, since it vanishes 
as k — > 0. We will show below that one must modify the spectral density for particles 
interacting with soft (non-hard-core) pair potentials because particle-shape information is 
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required in order to ascertain whether the system is hyperuniform or nearly hyperuniform. 
For particles interacting with a hard-core repulsion, the particle shapes are obviously the 
hard cores, but for non-hard-core interactions, such as in the Kob- Andersen model studied in 
this paper, one must determine a self-consistent procedure to assign particle shapes to each 
point particle. In addition, for such soft binary mixtures, the standard direct correlation 
function c(r), applicable to monodisperse systems, must be generalized. 

In this section, we present two generalizations of c(r) for polydisperse systems: one that 
is based on the spectral density (Sec. Ill BT]) . and another that is based on the matrix version 
of the structure factor (Sec. Ill B 2\\ . 



A. Single-Component Ornstein-Zernike Equation 

For a statistically homogeneous and isotropic single component system, the Ornstein- 
Zernike (OZ) equation 12 - defines the direct correlation function c(r) in term of the total 
correlation function h(r) = (?2( r ) — 1 and the system point density p: 

h(r) = c(r) + ph{r) <g> c(r), (6) 

where <g> denotes a convolution. After taking the Fourier transform of Eq. ([6]) and introducing 
the structure factor S(k) = 1 + ph(k), we get 

h{k) = c{k)+ph{k)c{k), (7) 

which is equivalent to expression (j2J). 

As noted in the introduction, we see from relation (T5]) that for hyperuniform single- 
component systems, i.e., lim^o S(k) = 0, c(k) diverges toward — oo as k approaches 0. 



B. Generalization of the Ornstein-Zernike Equation for "Two-Phase" Decorations 

As indicated in the beginning of the section, we must obtain a modified version of the 
direct correlation function for binary mixtures in which the particles interact with non- 
hard-core pair potentials in order to detect hyper uniformity or near-hyperuniformity. This 
function must be defined to be as general as possible. In particular, it must be equivalent 
to the usual direct correlation function in the case of a single- component system. We shall 



7 



therefore start by determining what this modified function would be in the single- component 
case in order to provide insight for the more general case of multiple-component systems. 
This will be done by decorating the underlying point configuration with nonoverlapping 
spheres. We first describe the single-component case and then the mixture case. 

1. Single- Component Case 

Consider a configuration of N points within a large volume V in which the minimum 
pair separation is the distance R. Now let us decorate this configuration by circumscribing 
spheres of radius R around each of the points, leading to a configuration of N nonoverlapping 
spheres of radius R. In this case, the particle phase indicator X(x) in terms of the positions 
of the sphere centers r 1; r 2 , . . . , rjv is:— >2i 

N 

J(x) = 5^m(|x-r i |; J R), (8) 

i=i 

where m(r; R) is the single-inclusion indicator function given by 

{1 r < R 
(9) 
0, r > R. 

The two-point correlation function S^r) = (-f( x K( x + r )) f° r such a statistically homo- 
geneous and isotropic distribution of nonoverlapping spheres, equal to the probability of 
finding two points, separated by the distance r = |r|, anywhere in the region occupied by 
the spheres, has been shown to be given by the following sum of two terms:— 

(f) — pm{r) <g> m(r) + p 2 m{r) <g> (?2( r ) ® m(r) , (10) 

where p = limy^oo N/V is the number density and angular brackets denote an ensemble 
average. The quantity pm®m is the self- correlation term, which is equal to the probability 
of finding two points inside the same sphere, and p 2 m £g> g% <g> m is the two-body correlation, 
the probability of finding two points in two different spheres. The autocovariance function 
X(r) is: 

x(r) = S*2(r) — p 2 v\ = pm(r) ® m(r) + p 2 m(r) ® ^(r) <S> m(r) — p 2 v\ ) 

= pm(r) <S> m(r) + p 2 m(r) <S> h(r) (g) m(r), (11) 
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where d 

Vl (R) = J m(r-R)dr = ^0 I - ) (12) 

is the volume of a (^-dimensional sphere of radius R [vi(R) = 47nR 3 /3 for d = 3 and V\(R) = 
ttR 2 for d = 2] . Taking the Fourier transform of Eq. (iTTj) yields 

X(k) = pm 2 {k) + p 2 m 2 {k)h{k) = pm 2 {k)S{k), (13) 

where S(k) is the structure factor define in (1T1). One can see from this equation that if the 
decorated "two-phase" nonoverlapping sphere system is hyperuniform, both S(k) and x(k) 
go to zero as k — > (phase in this context does not refer to a thermodynamical phase, but 
to either the particle or the void phase). 

In order to manage the extension of the standard direct correlation function that cor- 
respond to the autocovariance function x{k), we present the following analysis. The self- 
correlation term in relation must be subtracted because in its present form x( r ) is n °t 
analogous to h(r). Thus, we introduce a modified autocovariance H(r) = S2 (r) — pm(r) ® 
m(r), given explicitly by: 

H{r) = p 2 m{r) <g> h(r) <g> m(r). (14) 
Taking the Fourier transform of Eq. (1T41) leads to: 

H(k) = p 2 m 2 {k)h{k) = x(k) - pm 2 (k). (15) 

We can now define a new direct correlation function C(r) using H(r): 

H(r) = C(r) + Q(r)®C(r)®H(r), (16) 

where Q(r) is a function which is to be chosen such that lim^o C(^) diverges for any 
hyperuniform system, for which x(k) — > as k — > 0. 

H(k) = C{k) + Q(k)C(k)H(k), (17) 

d(k) = jW-^ 2 (fc) (18) 

l + Q(k) (x(k) - pm 2 (k)) 

For limk-^o C(k) to diverge for hyperuniform systems, we require that the denominator of 
the right side of Eq. f lTS]) to be zero whenever x(k) = 0, leading to the requirement 



Inserting Eq. ([19]) into Eq. ( ITT)) gives the one- component decorated OZ equation: 

H(k) = C(k) + ^Sg^, (20) 

= P m\k) - £^®. (21) 

Relation (I2TI) holds for a decorated single-component system. The generalization of 
Eq. fl2Tj) for a multiple-component system can be obtained by noting that Q~ l (k) is equal 
to the self-correlation term. For example, for a two-component system of nonoverlapping 
spheres, the relations analogous to f[TUj) - f[2"Tj) are given by 

^ W = p A m\(k) + p B m%(ky ^ (22) 

m = m - p A mi { k) - PB M = c ( k) + -^e^_ y (23 ) 

c(k) = + PB r n | W - (^iW^MM)! , (24) 

where and are the number densities of species A and B, respectively, and m^(r) and 
mg(r) are the corresponding sphere indicator functions. 



Mixture Case 



Consider an M-component system, in which N a represents the number of particles of 
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we write the following OZ equation for 



species a, where a = A,B,... Following Ref. 
the mixture total correlation function h a p(r) and the direct correlation function c a p(r): 

M 

h af s (r) = c a(3 (r) + ^ Pi c ^ ® h iP ( r ) > (25) 

7=1 

where a, (3, and 7 represent the different components of the system. Note that c a p{r) 
is different from the "decorated" "two-phase" direct correlation function C(r) defined in 
Sec. Ill B 1[ Equation (I25p can be rewritten in matrix form: 

(r) <g> ^p 1 pph lP {r) 

7 

H(r) = C(r) + C(r) g> H(r), (26) 
where the components of the matrices H(r) and C(r) are given by 

= y/p a pph a p(T), (27) 

Ca^(r) = y/pTp~pc a p{v). (28) 
10 



Taking the Fourier transform of Eq. (126]) gives 



H(fc) = C(fc) + C(fc)H(fc) 
C(Jfe) = H(ifc) (l + H(Jfe) 



(29) 
(30) 

where I is the identity matrix. 

Equation fl30l can be simplified by introducing the M x M multiple-component structure 
factor matrix S(k), whose components are denoted as S a p(k): 



( 



SAA(k) S AB (k) 
S\ B {k) S BB (k) 



\ J 
I + H(fc), 



\ 
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/ 



1 + PA^AA^k) a/ P aP E}h ab(Jz) 
y/pApBh* AB (k) l + p B h BB (k) 



V 



\ 
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(31) 



where denotes the complex conjugate of S a p(k). Substitution of Eq. f l3T|) into Eq. f )30|) 

yields the following simpler expression for C(k): 



C{k) =I-S(A;)- 1 . 



(32) 



This last equation should be used carefully, since the S(k) matrix is rank-1 for a single 
realization of a system, and hence it cannot be inverted without first taking an ensemble 
average.— 

Equation ff24|) . valid for the "two-phase" decoration, and Eq. (!30j) may not look similar, 
but their similarities can be made apparent by rewriting x(k) and S a p(\t) in terms of the 
collective coordinates p a (k): 

_A f Q _ 

(33) 



p a (k) = ^e^, 

where k is the wave vector and N a is the number of particles of species a For a single 
configuration of a multiple-component system in a volume V, we get the structure factor 
matrix components to be given by 



(34) 



Since we never compute S a p(k = 0) directly, instead relying on the k — > limit, we can 
drop the Kronecker delta function in the following steps. For a two-component system, the 
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spectral density for the decorated system is 



X(r) = p A m A (v) ®m A (r) + p B m B (r) <g)m s (r) + (35) 
p 2 A m A (r) <g) ^(r) <g) m j4 (r) + p A p B m A (r) g> h AB (r) ® m B (r) + 
PAPBm B (r) g) h BA (r) <8> m A (r) + p|m B (r) <g> ft, BB (r) <g> m B (r), 

for which the Fourier transform is given by 

X(k) = pA^i(k) + p B m 2 B (k) + p 2 A m 2 A (k)h AA (k) + p A p B m j4 (k)m B (k)/i j 4 i j(k) + 



p j4 p J Bm yl (k)m j B(k)/i j BA(k) + p B m B (k)h BB (k) 
\p A (k)m A (k) + p B (k)m B (k)\ 2 



V 

Using Eq. (J55|l to rewrite Eq. (J23|) leads to 



(36) 



G(k) = p x mi(k) + p B m|(k) f 1 Mm±Isnm \ 

V |p A (k)m A (k)+p B (k)m B (k)|V 

Now, assume that the decoration of the two-component system is chosen such that ^(k) = 
(y/pjfh^k), y/p^m B (k)) T is an eigenvector of S(k). Calculating the associated eigenvalue 
of C(k) (which shares eigenvectors with S(k)) leads to 

^* T (k)C(k)^(k) 1 N A m 2 A (k) + N B m 2 B (k) 



(38) 



p A m\(k) + p B m 2 B (k) (\p A (k)m A (k) + p B (k)m B (k)\ 2 ) ' 

The similarities between Eqs. fl37l) and fl38l) are striking, and lend credibility to their use. 
However, it should not be forgotten that Eq. (|38|) is only valid for a very precise choice 
of (k) and m B (k), which may or may not be realizable for arbitrary systems. It is 
therefore more appropriate to use a decoration that uses a priori information about the 
system (e.g. an effective radius of the particles) together with Eq. fl3Tl) . In a situation where 
such information is missing, calculating the actual eigenvalues of S(k) and C(k) is a good 
alternative choice, although it requires multiple realizations of the system in order to get 
the ensemble-average values. 



III. SIMULATION DETAILS 



We carry out molecular dynamics simulations in the NVT ensemble to study the behavior 
of two different atomic glass-forming liquid models: a three-dimensional single-component 
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system in which the particles interact with the Z2 Dzugutov potential and a two-dimensional 
two-component system in which the particles interact with the Kob-Andersen potential. In 
particular, starting from liquid states, we quench these two model systems and follow their 
transitions from fluids, to supercooled fluids and glassy states and determine their associated 
supercooled and glassy states as a function of temperature. 

The interacting systems consist of N = 100000 particles in a two-dimensional (Kob- 
Andersen) or three-dimensional (Z2 Dzugutov) periodic box, subject to a Nose-Hoover 
thermostat^ with a mass set to iV/1000 = 100. This particular choice of mass is selected 
to avoid the numerical instabilities that occur when a small mass is used, while reducing 
the time the thermostat takes to equilibrate which increases with larger masses. The initial 
configurations are generated using the random sequential addition (RSA) algorithm, 25 and 
with an initial temperature that is much larger than the freezing temperature. There are 
four relevant units in the molecular dynamics simulations: units of energy, length, mass, and 
time, of which three can chosen independently. The units of energy and length are selected 
by the numerical values of the potentials' parameters, while the unit of mass is set by letting 
all particles have unit masses. These choices defined the natural units, including the unit of 
time. The system is then continuously cooled using an exponential rate 

T(t) = T x 10~ t/ri °, (39) 

where T(t) is the temperature when the simulation has been running for a time t, T is the 
initial temperature, and the time per decade Tiq controls the cooling rate. The molecular 
dynamics integration is done using the velocity Verlet scheme. 

For the Z2 Dzugutov potential, shown in Eq. (|4]), we use the following parameter values: 
a = 1.04, rj = 0.33, k F = 4.139, b = 4.2 x 10 7 , a = 0.348, n = 14.5, r c = 2.64488, 
and Vq = 0.13391543. The values of r c and Vq are chosen such that both v(r c ) = and 
0. This choice of parameters defines the natural units of both energy and length. 



ill! 

dr 



r=r n 



Following Ref. |14j, the particle density is fixed as p = 0.84 and the particle mass is set 
to unity. The time per decade r 10 is set to 500, 200, and 50 natural time units. Slower 
cooling schedules are attempted (such as r w = 2000), but they lead to some of the samples 
crystallizing. The time step is At = 5 x 10~ 3 in the natural time units and is chosen such 
that the total energy of the system is conserved when the thermostat is removed. 

For the Kob-Andersen potential, shown in Eq. fl5]), we use a composition of particles 
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with number ratio A : B = 65 : 35 and the following parameters: Oaa = 1-0, €aa = 1-0, 
o~ab = &ba = 0.8, €ab = £ba = 1-5, (Tbb = 0.88, and = 0.5. The values for the Vo a p 
are chosen such that the potentials are continuous at r = 2.ba a p cutoffs. These choices of 
parameters define the natural units of energy (€aa) and length (ctaa)- Both particle species 
are assumed to have masses equal to unity. Following Ref. |2|, we set p = 1.161662. The 
time per decade of temperature decay r w is set to 2000, 400, 100, and 20. The time step is 
At = 1 x lO" 3 . 



IV. RESULTS 

A. Z2 Dzugutov Single-Component Glass 

To estimate the glass transition temperature T g of the Z2 Dzugutov model, we use the 
temperature at which the total energy per particle as a function of temperature changes 
slope most rapidly. Since the harmonic contribution to the average total energy per 

particle u has a constant slope, we subtract it from u to detect any change of slope. As 
seen in Fig. [TJ we obtain ksT g ~ 0.88 for the Z2 Dzugutov model. Comparatively, by 
observing the highest temperature at which the supercooled systems crystallized and the 
temperature at which such crystals melt, we roughly estimate the melting temperature to 
be T m /T g ~ 2.5 ± 0.5. 

To calculate the volume integral of the direct correlation function c(r), we need to find 
the limit of S(k) for k — > 0, and then substitute it in Eq. fl2]). Since S(k = 0) cannot be 
calculated directly in a finite simulation box of side length L because the smallest possible 
wavenumber accessible is 2n/L, an extrapolation from the available data to zero wavenumber 



must be used. Figure 2(a) shows the small-wavelength behavior of S(k) for Z2 Dzugutov 
model at different temperatures. It is clear that S(k) is nearly linear in k for k ^ 1, leading 
to a very good fit to a linear function. This linear behavior of S(k) for small k > implies 
that the real-space total correlation function h(r) decays, for large but finite r, as a power 
law — 1/r 4 or, equivalently, the direct correlation function decays as c(r) ~ — 1/r 2 . The 
numerical value of S(k = 0) only changes by up to 5% between the cubic fit for k < 2 shown 



in Fig. 2(a) and a linear fit for k < 1. Since the linear fit is less susceptible to overfitting 



and complex behavior for 1 < k < 2, we elect to use this linear fit to extrapolate the value 
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FIG. 1: (Color online) Strictly anharmonic portion of the total average energy (kinetic and po- 
tential) per particle u — 2>ksT of the system in term of the thermostat temperature T. This 
is obtained by averaging over 10 cooling simulations of supercooled Z2 Dzugutov systems using 
t\q = 400. 3ksT has been subtracted from the energy to help identify the glass transition. The 
glass transition temperature ksT g ~ 0.88 is estimated by finding the temperature at which the 
function slope changes most rapidly. The vertical dashed line is located at T = T g . The energy 
scale is normalized through our choice of potential parameters (see Sec. IIII|) , 
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(a) (b) 
FIG. 2: (Color online) Structure factors S(k) for Z2 Dzugutov systems supercooled using no = 500 



for various temperatures. The curves have been averaged over 10 realizations, (a) Cubic fits of the 
small- wavenumber (k < 2) structure factors. The type of fits and their cutoff are chosen such that 
they accurately reproduce the features of the structure factors, in particular the positive linear 
dependence near k = 0. [(b)] Larger- wavenumber structure factors. 
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of S(k = 0) for these systems. 
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FIG. 3: (Color online) Growing length scales for Z2 Dzugutov systems generated using various 
cooling schedules. For each cooling schedule, the results have been averaged over 10 realizations 
and fitted to the sum of an exponential and a linear function to smooth out the numerical noise. 



(a) Limit of c(k) for k — > 0, calculated using linear fits of S(k). |(b)| The static length scale £ c , 
defined by relation (|40|) . associated with these systems. Note that the nearest neighbor distance 
between particles at T = is 1.0539. 

From the Fourier transform of the direct correlation function c(k), which has units of 
volume, we define the following length scale: 

-c(0)] 1/d , 



(40) 



where d is the Euclidean dimension. From Fig. 3(a) , there is a striking evidence that c(k = 0) 
grows to a large negative value in the supercooled regime, leading to a doubling in the value 
of the length scale £ c . 

In the case of a single-component system at equilibrium, the compressibility relation links 

to its structure factor as follows: 



its isothermal compressibility k t = —y ^ 



purksT = S(0). 



(41) 



However, supercooled liquids and glasses are not equilibrium states and consequently 
Eq. ( )4T1) tends not to be satisfied. Following Ref. Q, we use the deviation from Eq. (HIT) 
to measure a nonequilibrium index X: 

S(0) 



X 



pn T k B T 
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1. 



(42) 




FIG. 4: (Color online) Nonequilibrium index X for Z2 Dzugutov systems supercooled using various 
cooling schedules denned in Eq. (jUJ). 



The isothermal compressibility k,t is computed by the following finite difference formula: 

AV 1 



V AP' 



(43) 



where AV^ is the change in volume of the simulation box and AP is the resulting change 
in pressure of the system after it is allowed to relax at constant temperature. The pressure 
is calculated using the virial relation. It bears mentioning that since the system is not at 
equilibrium, it is not in a steady state even before the change in volume. To minimize 
the impact of the uncompressed system relaxation, both the uncompressed and compressed 
systems are allowed to relax for the same amount of time before measuring their pressures. 

In the case of the Z2 Dzugutov system, we use a change of volume AV jV = 0.3% and 
the pressure is sampled from t — 5 to t — 10, where t = denotes the time at which the 
system is compressed. As can be seen in Fig. HJ X is zero^ for T/T g > 2, with only slight 
deviations due to noise and numerical inaccuracies. However, as the temperature is lowered 
to values approaching the glass transition, X increases up to a value of ~ 0.2 at T/T g = 1. 
For T < T g , the inability of the system to relax in a time of the order of the cooling schedule 
time per decade t w results in nearly constant values of k t and S(0) which leads to the 
asymptotic behavior of X as T — > 0. 

Is the growing nonequilibrium index X, a purely static quantity, correlated with the grow- 
ing relaxation times as the temperature decreases during the supercooling process? Figure [5] 
shows a positive correlation between X and r, where r is the timescale associated with the 
early relaxation process, extracted from an exponential fit function ~ e~ l l T of the system to- 
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FIG. 5: (Color online) Timescale r of the early relaxation process of the system versus the nonequi- 
librium index X. Both quantities have been averaged over 10 configurations. The circle are cen- 
tered on the averages of X and r, while the horizontal and vertical lines represent their respective 
uncertainties, with their half-length set equal to the average standard deviations. The initial con- 
figurations which are allowed to relax at constant temperature are generated from the liquid phase 
through a cooling schedule employing no = 50. Each datum represents a single temperature. 
Observe that r and X are positively correlated. Therefore, since X is a monotonically decreasing 
function of the temperature T (see Fig 2]), r also increases with decreasing T. The values of T/T g 
associated with each datum are, in order of smallest to largest r are as follows: 1.80, 1.61, 1.43, 
1.28, 1.14, 1.01, and 0.90. 

tal energy. To observe this process, we start with configurations that have been supercooled 
to a given temperature following a specific cooling schedule. These configurations are then 
allowed to evolve at constant temperature. It can be clearly seen that X and r are strongly 
and positively correlated. 

B. Kob-Andersen Aq^B-^ Two-Component Glass 

To calculate the spectral density x(fc), we decorate the systems by circumscribing disks 
of radius Ra and Rb centered around the point particles of species A and B, respectively. 
Since our derivation in Sec. Ill Bl requires the disks to be nonover lapping, we chose the largest 
possible radii that satisfy this condition. In the case of a Kob-Andersen glass, A particles 
are often located next to one another, while B particles can be further apart. This leads 
to our decision to use the distances between the closest A- A and A— B pairs of particles to 
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FIG. 6: (Color online) Example of a decorated Kob- Andersen glass configuration (a small subregion 
of the configuration only). The larger disks represent the A particles, while the smaller disks 
represent B particles. The radii of the disks are chosen such that the two closest A particles of 
the whole configuration are in contact and the closest A—B pair of particles are in contact. The 
configuration shown has been generated using no = 100, and is at a temperature of T/T g = 
6.7 x 10" 5 . The particle radii are R A = 0.513720 and R B = 0.329883 {R a /Rb = 1.55728). 



define the particle radii. Figure E] shows part of a glass configuration decorated using this 
procedure. 

In an identical fashion to the Z2 Dzugutov system, we use the change in slope of the total 
energy in terms of the temperature to estimate the glass transition temperature T g for the 
Kob- Andersen system. Since the Kob- Andersen system that we analyze is two-dimensional, 
its harmonic contribution to the energy is 2/cgT, which we subtract from the total average 
energy per particle u to detect any change of slope. The result obtained from Fig. [7| is 
T g ~ 0.31, which is reasonably close to the previously-reported value of T g = 0.33. 27 

As in the case of Z2 Dzugutov systems, the spectral densities x(fc) for Kob-Andersen 
liquids, supercooled liquids, and glasses have nearly linear behavior for k ^ 1. It is thus 
possible to prescribe a linear fit to extrapolate the values of x(k = 0), which is required to 
calculate C(k — 0) using Eq. (j2ljh We again define a length scale based on the C(k = 0): 

1 l / d 



c 



-C(0) 



(44) 



where d is the Euclidean dimension. Figure 9(a) shows the large change in value of C(k = 0) 
as the Kob-Andersen liquids are supercooled, leading to the length scale £c to increase by 
a factor larger than 5 between the fluid states and the zero-temperature glassy states. 
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FIG. 7: (Color online) Strictly anharmonic portion of the total average energy (kinetic and poten- 
tial) per particle u — 2UbT of the system in terms of the thermostat temperature T. This is obtained 
by averaging 10 cooling simulations of supercooled Kob- Andersen systems using no = 400. 2fc#T 
has been subtracted from the energy to help identify the glass transition. The glass transition tem- 
perature T g ~ 0.31 is estimated by finding the temperature at which the function slope changes 
the most rapidly. The vertical dashed line is located at T = T g . The energy scale is normalized 
through our choice of potential parameters (see Sec. IIIip . 
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FIG. 8: (Color online) Spectral density x(&0 versus wavenumber k for Kob-Andersen ^5-635 
systems supercooled using no = 400. The curves have been averaged over 10 realizations and 
fitted using fourth degree polynomials. The type of fits have been chosen for their ability to 
reproduce accurately the features of the structure factors for the range presented (0 < k < 3) . The 
disk radii for the decorations are calculated independently for each configuration. 
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(a) (b) 

FIG. 9: (Color online) Growing length scales for two-dimensional Kob- Andersen systems. For each 
cooling schedule, the results have been averaged over 10 realizations and fitted to the sum of an 



exponential and a quadratic functions to smooth out the numerical noise, (a) Limit of C(k) for 
k — > 0, calculated using the linear fits of |(b)| The static length scale £c, defined by relation 
|, associated with these systems. 




FIG. 10: (Color online) Smallest eigenvalue of lim^o C(k), calculated using a linear fit of the 
matrix structure factor S(/c). While the qualitative behavior of this eigenvalue can be compared to 



C(k = 0) (see Fig. 9(a) ), their quantitative values cannot directly be compared because they have 



different units: C(k) has units of volume, while C(k) is dimensionless. 



As mentioned in Sec. Ill B 21 there is a second generalization of the direct correlation 
function which does not require any a priori knowledge or about the particle shapes. Instead, 
one can use the matrix direct correlation function C(r) and its Fourier transform C(k). As 
can be observed in Fig. [TOj the qualitative behavior of the smallest eigenvalue of C(k) in 
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the k — > limit is strikingly close to the behavior of C(k) in the same limit. This indicates 
that our decoration choice is appropriate for detecting long-range density fluctuations in 
Kob- Andersen glasses and supercooled liquids. 
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FIG. 11: (Color online) Nonequilibrium index X for Kob- Andersen systems supercooled using 
various cooling schedules defined in Eq. (|47|) . 



Since the compressibility relation (|4ip applies only to single- component systems, we must 
generalize the nonequilibrium index X for mixtures. The compressibility relation for multi- 
component systems at equilibrium, is given by:— 



KTksT 



IBI 



M 
0=1 



IB 



(45) 



where the components B a » of the matrix B are 



B a n = — lim S a p(k), 



(46) 



|B| is the determinant of B, and |B| Wi g is the a (3 minor of B. The nonequilibrium index 
X for multicomponent systems can now be defined by using the mismatch between the left 
and right sides of Eq. ( l4"5j) . that is, 



X 



IBI 



K T k B T}2 a= l Efl=ll B 



M 



1. 



(47) 



a/3 



As for single-component systems, the isothermal compressibility for this multicomponent 
system is obtained by computing the virial pressure response to an incremental change in 
volume using Eq. (]43l) . 

For the Kob- Andersen system, we use a change of volume AV/V = 0.2% and the pressure 
is sampled from t = 20 to t = 40, where t — denotes the time at which the quenching is 
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halted and the system is compressed. As can be seen in Fig. [TT], X is zero^ for T > 2T g . 
Similarly to the phenomenon observed in the case of the Z2 Dzugutov system (see Fig. Sj), 
X increases up to a value of a value of ~ 0.15 at T = T g . The asymptotic behavior of X 
for T < T g is again the consequence of the system inability to relax in a time comparable to 
the cooling schedule time per decade r 10 . 

V. CONCLUSIONS AND DISCUSSION 

We have demonstrated here that the static structural length scales £ c and £c are able 
to distinguish subtle structural differences between glassy and liquid states, which extends 
the analogous results for metastable hard spheres^ to atomic thermal systems. Since these 
length scales are based on the volume integral of the direct correlation function c(r) and its 
generalization C(r), respectively, their growth as a liquid is cooled past its glass transition 
is a sign of the presence of long-range correlations in the glassy state that are not present in 
liquids. Additionally, the continuing increase of £ c and £c past the glass transition indicates 
that, while particles primarily undergo sequences of local rearrangements, the glass may 
still exhibit order on a significantly larger length scale as the system continues to cool. 
Our results using two-dimensional Kob- Andersen binary mixtures and three-dimensional Z2 
Dzugutov single-component systems, as well as the previous results for MRJ packings as 
evidence, we postulate that these length scales are relevant in various glasses. This includes 
not only atomic systems possessing pair potentials with steep repulsions and short-range 
attractions, but network glasses as well. For example, in a recent computational study,— 
which is supported by recent experimental results,= it was shown that realistic models of 
amorphous silicon can be constructed to be nearly hyperuniform, which implies that such 
glassy tetrahedrally-coordinated networks are characterized by a large static length scale 
£ c We also have shown that the nonequilibrium index X is positively correlated with a 
characteristic relaxation time scale, since they both increase as a system is supercooled. 

An interesting issue concerns the explication of the underlying geometrical reasons for 
the negative algebraic tail in the pair correlation function,— which also has been observed in 
hard-sphere systems.- (The former is exhibited for large but bounded pair distances, while 
the latter is valid asymptotically as r — > oo.) The local geometric diversity of particle 
arrangements in an amorphous solid medium inevitably creates short-range density fluctu- 
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ations. In particular, this is true for the nearly hyperuniform cases examined in this study. 
Without being too specific, one can formally divide a "jammed" particle configuration into 
two equal subsets containing particles experiencing either lower or higher local densities 
than the overall system average. The fact that the pair correlation functions display nega- 
tive algebraic tails with increasing separation r has basic implications for the relative spatial 
distributions of these low and high local density particles. In particular, it indicates that 
large numbers of either particle type cannot fit together to form arbitrarily large clusters 
that dominantly exclude the other particle types. Instead, their spatial patterns evidently 
involve interpenetrating percolating networks in three dimensions and highly non-convex 
clusters in two dimensions. The detailed statistical geometric description of these patterns 
and why they generate algebraic pair correlation function tails constitutes an important area 
for future investigation. 

The quantity X introduced earlier in Eq. (142]) as a measure of deviation from thermal 
equilibrium can be usefully interpreted in terms of system occupancy on the many-body po- 
tential energy landscape.— 1 Specifically, this focuses on the comparative behaviors of isother- 
mal compressibility at high-temperature thermal equilibrium in the liquid phase as opposed 
to the measured isothermal compressibility in the non-equilibrium glass phase in the T — > 
limit. In the former case, an incremental pressure change and accompanying volume change 
will include shifts in occupancy probabilities for the separate basins that tile the landscape; 
these shifts involve interbasin local particle rearrangements that act to enhance the volume 
change induced by the pressure perturbation. In contrast, at very low temperatures, the sys- 
tem is trapped in its initial basin; intrabasin vibrational motions have insufficient amplitude 
to allow the system to take advantage of the previous kinds of local particle rearrangements. 
The resulting absence of enhanced volume change due to those interbasin transitions reduces 
isothermal compressibility, causing X to increase above zero. 
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